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1  Abstract  or  Summary 

This  study  aims  at  developing  and  testing  anomalous  (turbulent)  transport 
coefficient  models  that  include  the  effects  of  plasma  turbulence  and  that  can 
be  used  in  more  realistic  numerical  simulations  of  MPD  thruster  flows.  The 
anomalous  coefficients  were  derived  using  a  kinetic  description  of  microinstabil¬ 
ities  in  the  finite-beta,  magnetized,  collisonal  and  flowing  plasma  of  the  MPD 
thruster.  Nonlinear  plasma  theory  was  lased  under  the  valid  approximation  of 
weak  turbulence.  The  resulting  transport  models  were  studied  in  detail  in  the 
parameter-space  corresponding  to  typical  MPD  thruster  operation  and  a  strong 
dependence  of  anomalous  transport  on  the  electron  Hall  parameter  was  found. 
The  models  were  than  reduced  to  a  set  of  compact  polynomials  ideally  suited  for 
self-consitent  fluid  flow  calculations.  The  usefulness  and  value  of  these  models 
were  demonstrated  by  including  them  in  an  advanced  plasma  fluid  code  and  car¬ 
rying  numerical  simulations  of  realistic  MPD  thruster  flows.  The  results  show 
the  impact  of  turbulence  on  various  aspectgof  dissipation  in  the  thruster  as  well 
as  on  the  overall  propulsive  efficiency. 
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2  Background 

2.1  Project  and  Final  Report  Description 

In  this  final  report  we  present  the  final  results  of  a  two-year  theoretical  and 
numerical  study  in  which  we  developed  and  tested  transport  coefficient  models 
that  include  the  effects  of  plasma  turbulence.  The  resulting  anomalous  transport 
models  are  intended  to  be  included  in  fluid  codes  thus  rendering  more  realistic 
the  numerical  simulations  of  MPD  thruster*  flows. 

Computer  models  of  the  MPD  thruster  rely  on  the  numerical  solution  of  the 
magnetohydrodynamics  (MHD)  equations  to  predict  the  performance  and  aid 
in  the  design  of  higher  efficiency  thrusters.  Current  MHD  models  of  the  accel¬ 
erator  include  only  classical  transport^  (i.e.  transport  due  to  collisions  between 
particles)  and  overpredict  the  performance  of  the  experimental  prototypes.  The 
inclusion  of  anomalous  transport  (i.e.  transport  due  to  “collisions”  of  particles 
with  the  oscillating  fields  of  unstable  waves  in  the  plasma)  is  expected  to  greatly 
enhance  the  validity  and  applicability  of  such  models.  During  the  past  year  our 
reserach  has  produced  a  first  generation  of  transport  coefficient  models  that  can 
be  readily  included  in  any  plasma  fluid  code  computer  code. 

2.2  Report  Organization 

This  final  report  documents  the  entire  work  starting  with  its  relevance  to  the 
design  of  high  efficiency  plasma  thrusters  for  spacecraft  propulsion,  the  method¬ 
ology  behind  the  approach,  the  derivation  of  the  anomalous  transport  models,  a 
description  of  how  they  are  to  be  used  in  a  standard  fluid  flow  code  and  ending 
with  actual  numerical  simulations  of  MPD  thrsuter  flows  using  the  models  in 
an  advanced  plasma  fluid  code  developed  at  our  laboratory. 

Non-technical  information  related  to  the  grant  and  the  work  it  supported  is 
given  at  the  end  of  this  report. 


2.3  Relevance 

The  ultimate  goal  of  MPD  thruster  research  is  the  development  of  high-efficiency 
long  lifetime  plasma  thrusters  that  can  be  used  for  energetic  spacecraft  propul¬ 
sion  requirements^. 

The  current  performance  of  MPD  thrusters  is  limited  to  modest  efficiencies 
(below  30%  for  argon)  and  plasma  fluid  computer  codes  (which  are  the  ultimate 
non-experimetal  tools  for  studying  thruster  design  problems)  are  handicaped 

'Throughout  this  report  the  term  MPD  thruster  refers  to  the  coaxial,  sdf-6ekl,  high-current, 
gas-fed  MagnetoPlasmaDynamic  thruster. 

^Throughout  this  report,  “transport"  refers  to  the  following  processes;  electrical  conductivity, 
particle  diffusion,  viscosity  and  thermal  conductivity. 

*By  “energetic"  we  mean  missions  that  require  high  increment  in  spacecraft  velocity,  Ae  (delta 
vee)  such  as  orbit  raising,  some  demanding  attitude  control  situations  and  interplanetary  propulsion. 
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by  the  absence  of  the  representation  of  some  crucial  physics.  Notably,  all  codes 
in  use  today  cannot  account  for  the  effects  of  plasma  turbulence  caused  by 
microintsbailies  in  the  plasma.  The  reason  such  turbulent  effects  are  important 
is  that  can  strongly  enhance  heating  and  ionization  of  the  plasma.  Heating  and 
ionization  are  essentially  dissipative  processes  for  the  MPD  thruster  because 
the  flow  is  practically  frozen  and  the  energy  invested  in  them  is  very  difficult  to 
recuperate  as  propulsive  thrust. 

Our  work  previous  to  this  AFOSR-sponsored  projectjl,  2,  3,  4]  has  shown 
that  the  MPD  thruster  is  a  plasma  accelerator  whose  performance  is  dictated 
by  the  competition  between  current-driven  acceleration,  which  constitutes  the 
prime  acceleration  mechanism  in  the  device,  and  current-driven  dissipation. 
Current-driven  acceleration  makes  the  MPD  thruster  a  prime  candidate  for 
spacecraft  propulsion  while,  current-driven  dissipation  is  a  handicap  that  must 
be  surmounted  if  these  devices  are  to  be  of  practical  use.  This  competition 
between  two  processes  driven  by  the  same  source,  the  current,  is  manifested  in 
the  low  thrust  efficiency  of  the  present  devices. 

The  high-current  acceleration  of  a  plasma  can  be  explained,  on  the  kinetic 
level,  as  primarily  due  to  the  transfer  of  momentum  from  the  current-carrying 
electrons  to  the  ion  gas  through  long-range  (Coulomb)  collisions.  Alternatively, 
on  the  macroscopic  level,  the  bulk  plasma  acceleration  can  be  viewed  as  the 
result  of  the  action  of  the  Lorentz  body  force  on  the  bulk  plasma.  The  Lorentz 
force  is  due  to  the  interaction  of  the  current  carried  by  the  electrons  with  the 
induced  magnetic  field. 

On  the  other  hand,  the  same  current  is  responsible  for  enhanced  (turbulent) 
transport  associated  with  current-driven  plasma  instabilities.  As  is  well  known 
from  plasma  physics,  any  current  having  an  associated  drift  velocity  exceeding 
the  thermal  velocity  of  one  of  the  charged  species  represents  a  potential  free 
source  of  energy  that  can  effectively  excite  one  or  more  of  the  many  natural 
plasma  wave  modes  into  unstable  growth.  For  the  MPD  accelerator,  the  nature 
of  the  current  (its  transversality  to  the  magnetic  field),  the  magnitude  of  its 
associated  drift  velocity  compared  to  the  ion  thermal  velocity  and  the 
plasma  parameters  are  such  that  a  current-driven  instability  (related  to  a  finite- 
beta  collisional  lower  hybrid  oscillations)  is  naturally  excited. 

This  current-driven  plasma  instability  has  the  peculiarity  of  anisotropic  elec¬ 
tron  heating  which  means  that  fast  electrons  are  preferentially  heated  along 
the  magnetic  field  lines.  Such  production  of  suprathermal  electrons  by  the  in¬ 
stability  usually  translates  into  a  substantial  enhancement  in  ionization  (over¬ 
ionization)  since  the  ionization  process  in  the  MPD  thruster  is  primarily  due  to 
electron  impact  i.e.  collisions  between  energetic  electrons  and  neutrals.  Sim¬ 
ilarly,  the  populations  of  excited  states  are  enhanced  (excitation)  through  the 
same  process. 

Plasma  heating,  excitation  and  over-ionization  are  not  desirable  in  the  MPD 
accelerator  since  the  time  scales  for  enthalpy  recovery  (  conditioned  by  the 
reverse  processes:  cooling,  deexcitation  and  recombination)  are  much  longer 
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than  the  plasma  residence  time.  The  translational,  internal  and  chemical  modes 
are  in  a  strongly  non-equilibrium  state,  or  equivalently  the  flow  is  frozen  with 
respect  to  these  modes  primarily  because  of  the  low  particle  densities  that  are 
associated  with  the  operation  of  MPD  devices  at  reasonable  power  levels.  In 
other  words,  the  low  density  typical  of  MPD  thrusters  limits  the  collisional 
equilibration  and  recuperation  of  the  energy  tied  in  the  random  kinetic,  internal 
and  chemical  modes  which  become  unavailable  to  the  propulsive  action  of  the 
plasma.  Thruster  efficiency  is  thus  degraded. 

The  work  documented  here  aims  at  developing  models  for  the  anomalous 
transport  caused  by  plasma  microinstabilities  as  well  as  at  the  testing  of  these 
models  in  actual  plasma  fluid  codes.  We  start  with  the  derivation  of  these 
models. 


3  Anomalous  Transport  Models  for  the  MPD 
Thrsuter  Plasma 

3.1  Methodolgy 

We  start  by  using  the  linear  stability  description  that  we  have  developed  in  work 
prior  to  this  AFOSR-sponsored  project[3, 4)  along  with  plasma  weak  turbulence 
theory  to  develop  a  second  order  description  of  wave-particle  transport  and 
anomalous  dissipation. 

In  Sections  3.2  and  3.3  we  outline  the  basic  formalism  we  adopt  for  our 
formulation  of  anomalous  transport.  In  Section  3.4  we  use  the  statistical  de¬ 
scription  of  the  previous  section  to  derive  general  finite-beta  expressions  for 
the  anomalous  ion  and  electron  heating  rates  as  well  as  for  the  electron-wave 
momentum  exchange  rate  that  controls  the  anomalous  resistivity  effect.  These 
expressions  are  cast  as  integrals  in  wavevector  space  of  quantities  that  depend 
on  the  various  elements  of  the  linear  dispersion  tensor  derived  in  detail  in  [5],  on 
the  roots  of  the  linear  dispersion  relation  and  on  the  saturation  energy  density 
of  the  fluctuating  (turbulent)  fields  denoted  by  €f 

We  then  turn  our  attention  in  Section  3.5  to  the  difficult  question  of  the 
saturation  mechanism  that  dictates  the  magnitude  and  dependencies  of  £f  We 
consider  four  models  for  £t  based  on  four  possible  saturation  mechanisms. 

In  Section  3.7  we  show  various  calculations  of  the  anomalous  heating  and 
momentum  exchange  rates  for  plasma  parameters  of  interest  and  compare  their 
magnitudes  to  classical  values. 

We  Anally  conclude  the  analysis  of  anomalous  transport  in  Section  4  by  using 
these  calculations  to  arrive  at  polynomial  expressions  of  the  relevant  transport 
coefficients  cast  solely  in  terms  of  macroscopic  parameters  for  inclusion  in  plasma 
fluid  codes  of  the  accelerator. 

These  polynomials  are  then  ready  for  use  in  a  fluid  code  as  described  in  the 
rest  of  the  report. 
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3.2  Weak  Turbulence  Kinetic  Formalism  for  Anomalous 
TVansport 


The  kinetic  theory  of  weak  turbulence  was  first  developed  by  Vedenov,  Velikhov 
and  5a9dcev[6]  (1961)  as  well  as  Drummond,  Pines  and  Rosenbluth  [7]  (1962). 
Its  current  status  is  probably  best  exposed  by  Galeev  and  Sagdeev  in  ref.  [8] 
(1982).  The  theory  allows  for  far  deeper  insight  into  the  description  of  turbulent 
transport  than  that  afforded  by  the  hydrodynamic  theory  of  turbulence.  This 
is  due  to  the  fact  that,  unlike  a  neutral  fluid  where  turbulence  is  synonymous 
with  strong  nonlinear  couplings,  plasma  microturbulence  is  often  diaracterized 
by  primarily  linear  interactions  between  unstable  modes  superimposed  on  small 
amplitude  nonlinearly  induced  waves[9,  p.  289). 

Formally,  the  validity  of  the  weak  turbulence  theory  hinges  on  the  smallness 
of  the  ratio  of  the  fluctuating  energy  density  to  the  plasma  thermal  energy 


density 


Si 


(1) 


where  n,  and  T,  are  the  density  and  temperature  of  species  s.  We  shall  find 
it  more  convenient  in  the  following  sections  to  deal  with  the  ratio  ft/noT*,  (r»o 
representing  the  background  charged  particle  density)  which  is  very  close  to  the 
above  ratio  for  a  quasi-equithermal  plasma.  It  is  possible  to  relate  the  weak 
turbulence  scaling  parameter  St/noTi  to  the  experimentally  measurable  density 
fluctuation  n/no,  where  the  tilde  denotes  a  fluctuating  quantity  by  noting  that 
n/no  w  c^/r«  and  ci>  «  eE/k  (where  e  is  the  electron  charge,  ^  the  wave 
potential,  k  the  wavenumber  or  the  magnitude  of  the  wavevector  k  and  E  the 
electric  field)  so  that 


Si  T,  jkr^f  /  n  \ 

noT.  ~  T,  4  Vno 


Experimental  evidence  of  turbulent  fluctuations  caused  by  cross-field  current- 
driven  instabilities  was  recently  found  in  the  low-power  steady-state  MPD  thruster 
plasma  at  various  conditions  and  locations  in  the  plume|10, 1 1].  These  measured 
turbulent  fluctuations  had  most  of  their  power  in  the  lower  hybrid  mode  with 
some  power  appearing  sometimes  in  the  electron  cyclotron  harmonics.  Mea¬ 
sured  values  of  (u/no)^  when  such  turbulence  was  observed  were  on  the  or¬ 
der  of  .1  with  magnitudes  ranging  between  .05  and  .7.  Fbr  these  values,  with 
1  >  Tj/T*  <  6,  ft  =  co|EtP/2  (where  cq  is  the  permittivity  of  free  space  and 
Tee  is  the  electron  cyclotron  radius)  and  assuming  {kra)"  w  .1(5],  we  obtain 
from  Eq.  (2)  an  estimate  for  £t/noTi  ranging  between  10~®  and  10“*  implying 
that  the  weak  turbulence  assumption  is  generally  valid. 

One  of  the  most  fortunate  turn  of  events  in  modern  plasma  physics  has  been 
the  surprising  ability  of  weak  turbulence  theory  to  describe  experimental  obser- 
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vations  even  for  fairly  large  amplitude  waves{9,  p.  291].  A  historical  perspective 
on  the  evolution  of  the  formalism  is  given  in  (5,  pp.  118-119). 


3.3  Governing  Equations 

3.3.1  The  Moment-Generating  Equation 


In  this  sub-section  we  present  an  outline  of  the  derivation  of  the  general  form  of 
fluid-like  equations  governing  the  evolution  of  macroscopic  quantities  under  the 
conditions  of  weak  turbulence.  Detailed  discussions  of  such  a  derivation  have 
already  been  presented  in  numerous  articles  (see  refs.  (12,  13],  for  instance,  for 
a  tutorial  review).  We  do  this  in  preparation  to  our  derivation  of  anomalous 
transport  presented  in  the  following  sections. 

We  should  mention  at  the  outset  that  our  interest  lies  not  in  the  evolution 
equation  itself  but  rather  in  its  use  as  a  moment-generating  equation.  Therefore, 
for  the  sake  of  simplicity  and  in  order  to  keep  a  connection  with  the  literature,  we 
shall  neglect  collisions  in  the  kinetic  evolution  equation.  The  effects  of  collisions 
will  be  reintroduced  later  when  we  use  the  explicit  form  of  the  dispersion  tensor 
elements. 

The  underlying  idea  is  to  consider  the  distribution  function  of  the  s  species, 
/,,  as  the  sum  of  a  slowly  varying  ensemble- average  part  and  a  rapidly  varying 
fluctuating  part 

f,{x,v,t)  =  F,iv,t)  + ft{x,v,t)  (3) 

where  F,(v,t)  =  {f,(x,v,l)}  and  {  >  denotes  an  ensemble-avBrage[14]  while 

the  tilde  denotes  a  quantity  fluctuating  due  to  the  effects  of  unstable  waves. 
When  similar  partitions  are  effected  on  the  electric  and  magnetic  field  vectors, 
the  kinetic  (Vlasov)  equation  for  a  spatially  uniform  equilibrium  yields 


(4) 


where  q,,  m„  iva  are  charge,  mass  and  cyclotron  frequency  of  species  ?  and 
where  we  have  chosen  to  work  with  the  cylindrical  phase  space  coordinates 
v±,4>,Vt  (representing  the  velocity  perpendicular  to  the  magnetic  field  B,  the 
cylindrical  angle  and  the  parallel  velocity,  the  magnetic  field  being  aligned  with 
the  2-axis).  Taking  the  ensemble-average  of  the  above  equation,  while  noting 
that  (/,)  =  0,  results  in 


dFs  dF,  _  (dF,\ 
dt  dl 


(5) 
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where  the  right  hand  side  represents  the  anomalous  contribution  that  is  the 
response  of  the  average  distribution  function  to  the  microturbulent  fluctuations 
and  can  be  written  explicitly  as 

By  subtracting  Eq.  (6)  from  Eq.  (4)  and,  in  the  spirit  of  weak  turbulence  the¬ 
ory,  neglecting  all  terms  that  are  quadratic  in  the  fluctuation  amplitude  (which 
is  tantamount  to  the  neglect  of  nonlinear  wave-particle  and  wave-wave  inter¬ 
actions)  the  following  governing  equation  is  obtained  for  a  weakly  turbulent 
plasma 

^  +  .  +  V.F,.  (T) 

The  standard  procedure  in  weak  turbulence  theory  (expounded  in  ref.  [14] 
for  instance)  is  to  solve  Eq.  (7)  along  with  Maxwell’s  equations  for  E,  B  and 
f,  then  substitute  the  result  into  Eq.  (5)  to  obtain  the  evolution  of  F,  in  the 
presence  of  microturbulence.  We  shall  not,  however,  need  to  do  all  that  for 
our  particular  problem  of  deriving  expressions  for  the  momentum  and  energy 
exchange  rates.  Such  expressions  can  be  arrived  at  by  taking  moments  of  the 
governing  equation  (Eq.  (5))  as  outlined  below. 


3.3.2  Evolution  of  Average  Macroscopic  Properties  under  Microtur- 
bulcncc 


To  obtain  the  macroscopic  evolution  equations  we  take  moments  of  Eq.  (5)  i.e. 
we  multiply  the  equation  by  the  generic  quantity  of  transport  0  (which  could 
represent  mass,  momentum  or  energy)  and  integrate  over  velocity  space  to  get 


dt 


/ 


QF,d  V  -  q,m. 


X  Bo  ■  V„0)  F,  dv  = 


TTl^ 


Q[(^E  +  vxBy\/„e\f,dv^ 


(8) 


where  we  have  used  integration  by  parts  in  order  to  move  the  distribution  func¬ 
tions  outside  the  operators.  Taking  successive  moments  of  Eq.  (5)  is  equivalent 
to  substituting  0  =  l,r,wu  (for  mass,  momentum  and  energy  respectively)  in 
Eq.  (8)  and  integrating  over  u-space.  This  yields 


d{n,} 

dt 

d{r.) 

dt 


=  0 

"h  (<*^C5®j)  ^  {E »)  = 


^(Eh,  +  r,xB) 
m,  \  / 


(9) 


(10) 
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^^  +  2(a;«e,)x(W,)  = 

2^(Er,  +  W,xB)  (11) 

where  c,  is  the  unit  vector  along  the  2-axis,  and  we  have  used  the  following 
definitions 


for  the  average  number  density,  the  particle  flux  density  and  the  kinetic  energy 
density  tensor,  respectively  (with  Vd,  as  the  drift  velocity  vector  of  species  s). 


3.4  Momentum  Exchange  and  Heating  Rates 

We  now  proceed  to  define  and  derive  explicit  relations  for  the  anomalous  rates 
of  interest. 

The  right  hand  side  of  Eq.  (10)  represents  the  rate  of  momentum  exchange 
{dPs/dl)AN  (where  the  momentum  density  vector  is  P,  =  m,r,)  between 
the  particles  and  the  fluctuating  fields.  Since  we  shall  be  interested  in  the 
momentum  exchange  along  the  drift  velocity  vector,  we  write 

=  (15) 


where  we  have  defined  as  the  effective  anomalous  momentum  exchange 

rate  (or  frequency)  between  species  s  and  the  fluctuating  fields.  Using  the 
explicit  expression  for  {dP,/dt)AH  from  Eq.  (10)  in  the  above  equation  we 
obtain 


Was  =  - 


n,Tn,Vds 


E  ■  w^n,  ^  Vdj  ■  {f,  X  D) 


Vd, 


(16) 


where,  unlike  most  derivations  in  the  literature,  we  are  retaining  the  full  elec¬ 
tromagnetic  character  of  the  microturbulence. 

We  now  specialize  the  above  expression  for  our  particular  problem  according 
to  the  MPD  thruster  configuration  shown  in  Fig.  (1).  We  thus  obtain  the 
effective  anomalous  momentum  exchange  rate  for  electrons  along  the  current 
after  setting  s  =  c,  staying  in  the  ion  reference  frame  and  aligning  the  relative 
drift  Ude  along  the  y-axis. 


Wan  = - - - +  noUde.Bz  -  ncUdt.B,)  (17) 

nomeU^e  '  ' 
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Figure  1;  The  vectors  j,  B,  k  and  u*.  in  the  local  cartesian  coordinate  frame. 
Also  shown  is  the  accelerator’s  fixed  cylindrical  coordinate  frame,  t-z'-O. 


where  we  have  used  the  relation  r,  =  +  noVi,. 

The  frequency  can  be  thought  of  as  an  effective  “collision”  frequency 

between  the  electrons  and  the  fluctuating  fields  and  can  thus  be  associated 
with  a  resistivity  called  “anomalous  resistivity”  the  same  way  that  the  coulomb 
collision  frequency  u^i  is  associated  with  the  classical  Spitzer  resistivity.  By 
analogy  the  anomalous  resistivity  is  proportional  to  ("Daw  and  is  given 


by 


{'1)as  = 


noe* 


(18) 


The  effective  collision  frequency  is  therefore  a  direct  measure  of  anoma¬ 

lous  resistivity. 

Similarly,  for  the  temperature 


T,  =  ^  J{v  -  Vitff.  dv 

we  define  a  heating  rate  for  species  s  by 


(19) 

(20) 


and  obtain,  after  combining  Eqs.  (10)  and  (11)  and  specializing  for  the  MPD 
thruster  configuration, 


3^  ("'■*  ■**) 


(21) 
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for  the  ions  and 


-2c 


-  /no£?  • 

e  ' 


tide 


ZnoTe 

noUje  («de.  Bx  -  tide.  B*)  ^ 


(22) 


for  the  electrons. 

Eqs.  (17),  (21)  and  (22)  will  be  the  focus  of  our  remaining  analysis  and 
calculations. 

In  order  to  proceed  with  more  useful  forms  for  these  expressions  we  need 
to  eliminate  the  fluctuating  density,  velocity  and  magnetic  field  in  favor  of  the 
fluctuating  electric  field.  For  this  purpose  we  invoke,  in  the  spirit  of  a  quasilin- 
ear  description,  relations  between  the  fluctuating  quantities  and  the  fluctuating 
electric  field  that  follow  those  of  their  linearly  oscillating  counterparts.  FVom  a 
generalized  Ohm’s  law  we  can  write  for  species  s[4) 

j,^  =  (23) 


(where  j  is  the  current  density,  is  the  dispersion  tensor  of  species  s  and  u> 
is  the  wave  frequency)  which,  combined  with  the  continuity  relation[4] 


k » 


(24) 


gives  a  useful  expression  for  the  fluctuating  density  of  species  s 


-  _  fc  j.  _ 

<^q5 


M  £• 

■Im 


(25) 


where  the  subscript  fc  is  a  reminder  that  these  relations  are  for  the  spectrally 
resolved  {i.e.  Fourier  transformed)  fluctuations.  In  this  expression,  are 
the  elements  of  the  tensor  representing  the  dielectric  response  of  species  s  and 
can  be  readily  obtained  from  Eqs.  (90)-(97)  (derived  in[5]  and  quoted  in  the 
appendix)  through  transformations  that  will  be  described  further  below. 

In  a  similar  fashion  we  can  derive  an  expression  for  ids  from  the  following 
relation 

3  s  =  9*  {rtQUds  +  ifUds)  (26) 

and  Eq.  (25),  yielding 


We  shall  not  need  to  worry  about  the  second  term  on  the  right  hand  side  of 
the  above  equation  in  the  context  of  the  MPD  thruster  configuration  shown  in 
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Fig.  (1),  because  this  term  vanishes  for  the  ions  =  0,  having  chosen  to  stay 
in  the  ion  rest  frame)  and  for  the  electrons  it  is  also  nil  for  the  components  that 
figure  in  Eqs.  (17)  and  (22)  (i.e.  the  x  and  z  components)  so  that  we  are  left 
with 


U)B 


(28) 


where  I  =  x, z  for  s  =  e;  and  I  =  x,y,z  far  s  =  i. 

Having  related  the  fluctuating  density  and  velocity  to  the  fluctuating  electric 
field  we  need  to  do  the  same  for  B.  To  this  end,  the  following  equation 


E  =  uajA  —  ik^. 


(29) 


which  relates  the  electromagnetic  potential  A  and  electrostatic  potential  to 
the  electric  field  E  (see  our  previous  study  in  (4)),  gives,  for  our  particular 
configuration. 


=  iujAr 


-'Vk 


K± 


£**  =  -ikt^k  + 


(30) 

(31) 

(32) 


Furthermore,  combining  the  above  equations  with  the  definition  of  the  electro¬ 
magnetic  potential, 

B  =  V  X  >4,  (33) 

and  Coulomb’s  gauge,  yields  the  desired  relations 


Bix  =  ^  (kxE^  - 


B. 


Vk 


B, 


-  —E 

=  -^4 


(34) 

(35) 

(36) 


We  are  now  in  a  position  to  evaluate  the  terms  of  Eqs.  (17),  (22)  and  (21)  by 
carrying  the  ensemble-averages  using  the  random  phase  approximation  (which 
is  a  standard  technique  of  statistical  physics  commonly  used  in  the  spectral 
resolution  of  fluauations,  see  Ref.  (15,  pages  371-373]  for  instance).  For  the 
first  term  in  Eq.  (17)  using  Eq.  (25)  we  have. 


>(e) 

Hm 


(37) 
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which  yields  under  the  assumption  of  random  phase^ 

(4n.)  =  -^j  dk  (38) 

where  9{  }  denotes  the  imaginary  part  of  a  complex  quantity.  Similarly,  we 

find  for  the  other  two  terms  in  Ekj.  (17) 

X  dk  (39) 

X  fcj.Ei|  dfc.  (40) 


If  we  now  substitute  the  above  three  equations  in  Eq.  (17),  expand  and  collect 
the  terms  in  the  summations  while  taking  adrantage  of  the  following  symmetry 
properties  of  the  dispersion  tensor 


we  arrive  at 


p(s)  -  _/?W-  ffW  =  -/?W.  oW  _  /?(•») 


*V 


V* 


«0 


UieTTlene 


jdfe. 


(41) 


(42) 


We  shall  find  it  convenient,  for  our  particular  instability,  to  cast  the  the  above 
expression  in  terms  of  the  spectrally  resolved  fluctuating  field  energy  density 
in  the  perpendicular  direction,  Sk^.  This  can  be  done  by  using  the  following 
relations  obtained  from  Eq.  (41)  and  the  governing  equation  [4] 

=  0  (43) 


(where  the  superscript  denotes  the  oscillating  part  of  the  electric  field)  to  yield 

A  =  ^  ~  RzyRyt 

“4  Rzytizx  +  RzzRyx  ^  ’ 


^FVoin  here  on,  we  shall,  for  the  sake  of  simplicity,  drop  the  subscript  k  from  the  fluctuating 
quantities. 
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D  —  ^  KryKxt  ~  RnRyx  , 

~R,.R»  +  RsMRtM 

(where  eadi  element  Rim  is  the  sum  of  the  rorresponding  contributions  from  the 
electrons,  ions  and  vacuum)  to  eliminate  and  E,  and  give 

UjeTTlene  J 

+Ei*>fl^  +  2fiiVB}dik.  (46) 

We  have  carried  out  our  derivation  above  under  the  electric  field  formalism 
where  the  relevant  dispersion  tensor  is  R  (cf  Eq.  (43)).  The  dispersion  tensor 
D  that  we  derived  explicitly  in  [5],  however,  was  obtained  under  the  potential 
formalism  (c/  Eqn.  (30)  of  ref.  [4]).  As  was  the  case  in  that  paper,  switching 
to  the  potential  formalism  has  some  advantages.  In  the  context  of  anomalous 
transport,  the  potential  formalism  allows  a  more  physical  insight  by  expressing 
the  momentum  exchange  and  heating  rates  in  terms  of  an  electrostatic  contri¬ 
bution  plus  a  finite-beta  correction.  The  results  obtained  so  far  can  readily  be 
recast  in  terms  of  the  elements  Dim  of  Eqs.  (90)-(97),  through  the  following 
linear  transformations  obtained  from  combining  Eqs.  (43),  (30),  (31)  and  (32), 


Ext  =  i?22  (47) 

=  -ft,*  =  y  (y£>23  -  A2)  (48) 

Rz,  =  -««  =  -^£>23  -  (49) 

Ryy  =  (50) 

+^{Du-Dj3)  (51) 

=  ^D^  +  !^Du-2^^Dn.  (52) 


We  also  need  to  separate  the  contributions  of  electrons,  ions  and  vacuum  in  the 
dispersion  tensor,  which  can  be  done  following 


/■)(«) _  n 

^Im  ~  ~  Dim  ~  D, 


where  subscripts  I  and  m  cover  the  indices  1,2  and  3  and  where  the  superscript 
(0)  denotes  the  contribution  of  vacuum.  The  elements  dI^  and  Dili  are  given 
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by 


and 


=  1. 

=  1  -  N® 

n(0)_ 

i/jj  — 

Df3>  =  Dg>  =  0 

(54) 

£)(*)  _ 

TiA 

p^(l+Ci^.) 

(55) 

II 

(56) 

£)(•)  _ 
Ui2  — 

o 

II 

II 

(57) 

When  the  above  transformations  (Ekjs.  (47)-(57))  are  used  in  Eq.  (46)  to  elimi- 
?im  favor  of 

2 


nate  in  favor  of  we  finally  obtain  after  some  straightforward  algebra 


~  li^e) AnIl  ~ 


uj^rrifrit 


+D«/l'  +  Og>[g  +  §e>-2S^ 


+2D'C^ 


B- 


1.2  t.2 

''J.*'*  r>2  , 

+ 


(58) 


where  [(t'J)y(^]i,  is  the  well-known  electrostatic  (longitudinal)  contribution  to 
the  anomalous  electron  momentum  exchange  rate 


M) 


AHiL 


Ujemefie 


f  dk 


(59) 


and  Xe  =  is  the  electrostatic  susceptibility  of  the  electrons.  In  the  electro¬ 
static  limit  (/3  — »  0)  it  can  be  verified  that  the  integrand  in  Eq.  (58)  vanishes  so 
that  we  are  left  with  (t'J)^,v  and 

Q{X«}  =  =  9{-x.}-  (60) 


We  shall  demonstrate  through  the  calculations  of  Section  3.7  that  the  trans¬ 
verse  (electromagnetic)  or  finite-beta  correction  to  ^ 

substantial,  especially  for  a  finite-beta  plasma  like  that  of  the  MPD  thruster. 

Equations  (46)  and  (58)  are  equivalent  but,  for  the  present  analytical  discus¬ 
sion,  we  prefer,  from  here  on  to  use  the  former  (ie.  the  electric  field  formalism) 
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because  the  resulting  expressions  are  more  compact.  For  our  numerical  calcu¬ 
lations  we  shall  apply  the  transformations  in  Eqs.  (47)  to  Eq.  (57)  in  order  to 
obtain  the  R  tensor  from  the  D  derived  in  our  previous  work  [4]  and  quoted 
here  in  the  appendix. 

It  is  convenient  to  express  (i^Dx/v  ’^™ts  of  a  natural  frequency.  We  choose, 
the  lower  hybrid  frequency,  wja,  because  it  is  the  natural  plasma  oscillation  mode 
closest  to  the  modes  we  are  investigating,  and  normalize  Eq.  (46)  to  get 

K'L'tJ  rrit  J  rtcTi 

^  2R«r]  }  d  fc.  (61) 

We  have  focused,  above,  on  the  anomalous  electron  momentum  exchange 
rate.  Similar  derivations,  with  no  conceptual  difference,  start  from  Eqs.  (21) 
and  (27)  and  lead  to  the  following  expressions  for  the  ion  and  electron  heating 
rates  and 


The  above  three  equations  are  the  sought  expressions  for  our  analysis  of 
anomalous  transport®. 

^Winske  et  o/.(16|  (1985)  have  derived  expressions  for  the  anomalous  heating  rates  of  ions  and 
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3.5  Saturation  Mechanisms 

For  the  numerical  analysis  of  anomalous  transport  in  the  MPD  thruster  plasma, 
the  above  three  equations,  along  with  the  tensor  elements  in  Eqs.  (90)-(97),  the 
linear  dispersion  relation 

detl/Zy(a;,fc)|  =  0  (64) 

and  the  transformations  in  Eqs.  (47)-(57)  form  an  atmoat  complete  set  of  equa¬ 
tions  in  terms  of  the  following  non-dimensional  parameters; 


U) 

Wifc’ 


7  Tt  Ude 

•  rr>  *  t 

Tg  Vti 


nij  Ve 

nit '  wjft 


(65) 


where  uJih  ~  y/UJ^iUZ  is  the  lower  hybrid  frequency,  0t  is  the  electron  beta  (ratio 
of  electron  thermal  pressure  to  magnetic  pressure),  Wpt  is  the  electron  plasma 
frequency  and  4*  is  related  to  the  propagation  angle  0  (see  Fig.  (1))  scaled  by 


the  mass  ratio 


(66) 


The  only  lacking  equation  is  one  that  relates  the  level  of  saturated  microturbu¬ 
lence  Sk/noTi  to  the  above  parameters. 

The  formulation  of  this  relation  is,  by  far,  the  most  formidable  obstacle  to 
the  study  of  wave-particle  transport,  since  such  a  relation  must  embody  all  the 
physics  of  the  nonlinear  saturation  mechanism.  The  central  question  in  this 
context  concerns  the  mechanism  through  which  the  fluctuations,  initiated  by 
the  instability,  reach  a  steady-state.  This  mechanism  dictates  the  magnitude 
and  dependencies  of  the  corresponding  fluctuating  energy  density. 

In  this  report  we  shall  adopt  and  compare  four  different  saturation  mecha¬ 
nisms;  ion  trapping,  electron  trapping,  ion  resonance  broadening  and  thermo¬ 
dynamic  bound.  The  models  for  each  of  these  mechanisms  are  discussed  in  more 
detail  in  Ref.  [5,  pp. 128-135].  Fbr  our  purposes  here  we  only  quote  the  resulting 
expression  for  each  of  these  saturation  models. 


electrons  in  finite-bets  plasmas.  Their  expressions  differ  from  ours  because  of  a  difierence  in  evalu¬ 
ating  the  ensemble-averages.  In  going  from  Eq.  (37)  to  Eq.  (38)  above,  for  instance,  these  authors 
would  neglect  in  the  summation  the  cross-terms  with  I  ft  m  (see  Equations  (14)  and  (IS) 
of  that  paper).  There  does  not  seem  to  be,  however,  a  valid  a  priori  reason  to  neglect  such  cross- 
terms[17].  Indeed,  had  we  dropped  these  terms,  the  last  term  of  the  sum  in  the  integrand  of  Eq.  (46) 
would  vanish.  This  term,  2R^}b^,  was  found  to  be  important  for  many  situations  we  have  investi¬ 
gated  numerically.  Furthermore,  had  we  dropped  the  cross-terms  it  would  not  have  been  possible 
to  recover  the  electrostatic  contribution  from  Eq.  (58)  in  the  /f  ->  0  limit. 
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3.6  Thermodynamic  Limit:  The  Fowler  Bound 

An  upper  limit  for  £t  was  first  derived  by  F<m>/er{18]  (1968)  from  thermodynamic 
arguments, 

(67) 

and  simply  states  that  the  energy  in  the  turbulent  fields  cannot  exceed  the 
kinetic  energy  of  the  electron  drift  that  is  fueling  the  instability.  Fbr  a  convenient 
incorporation  in  our  particular  formulation,  we  recast  the  inequality  above  to 
read 


3.6.1  Saturation  by  Ion  Trapping 

In  situations  where  the  excited  wave  spectrum  is  narrow  due  to  the  dominance 
of  a  single  wave  mode,  a  monochromatic  wave  saturation  model,  such  as  that 
behind  particle  trapping,  can  prove  to  be  a  viable  mechanism  for  saturation. 
In  such  a  case  the  saturation  dynamics  can  be  governed  by  the  trapping  of  the 
particles  in  the  potential  wells  of  the  growing  mode  thus  limiting  its  growth.  At 
saturation  one  can  simply  write 

=  irui  (y)  (69) 

where  Wr/k  is  the  phase  velocity  of  the  dominant  mode  and  we  have  assumed 
that  the  ions  are  the  particles  being  trapped.  Again,  we  normalize  the  saturation 
model  for  compatibility  with  transport  theory  so  that,  using  e<^  w  eE/k  and 
=  «o|5'fcP/2i  the  above  equation  can  be  rewiitten  as 


3.6.2  Saturation  by  Electron  IVapping 

Electron  trapping  is  probably  not  a  viable  saturation  mechanism  for  an  insta¬ 
bility  in  which  electrons  are  collisional  and  are  in  broad-resonance  with  the 
unstable  waves.  We  shall,  however,  include  a  model  for  electron  trapping  satu¬ 
ration  in  our  calculations  for  the  sake  of  reference.  In  analogy  with  ion  trapping, 
we  can  write  for  the  electrons  as  viewed  from  the  ion  rest  frame 


and  after  some  algebraic  manipulations, 
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3.6.3  Saturation  by  Resonance  Broadening 

This  mechanism  relies  on  the  broadening  of  wave-particle  resonances  by  the 
random  motion  of  particles  in  the  turbulent  electric  held  set-up  by  the  microin¬ 
stability. 

If  resonance  broadening  is  to  be  important  in  our  case,  it  would  most  prob¬ 
ably  rely  on  ion  dynamics,  since  the  electrons  are  already  broadly  resonant 
with  the  waves  due  to  collisions  and  finite-beta  effects  while  the  ions  are  very 
narrowly-resonant(5] . 

Following  Gary  and  Sandcr5on(19]  who  applied  the  Dum-Duprte  resonance 
broadening  formula[20]  to  the  ions  and  found,  after  taking  the  velocity  average 
/  Aw/i,  dvf  /  /i,  dv, 

(^fc);RB  =  2^°^°  (x)  ’ 

we  specialize  the  ion  resonance  broadening  model  for  our  dimensionless  param¬ 
eters  and  obtain 

/  u>r  y  .... 

noTi  (ifcr^f  Ti  W/  U/J  ■  ^  ’ 

3.7  Calculations  of  Anomalous  TVansport 

Armed  with  the  expressions  for  anomalous  transport  in  Eqs.  (61),  (62)  and 
(63)  along  with  the  tensor  elements  in  Eqs.  (90)-(97)  (Appendix),  the  linear 
dispersion  relation  in  Eq.  (64),  the  transformations  in  Eqs.  (47)-(57)  and  the 
saturation  models  in  Eqs.  (68),  (70),  (72),  (74),  we  can  now  conduct  a  compar¬ 
ative  numerical  study  of  anomalous  dissipation. 


3.7.1  Classical  Benchmarks 


For  benchmarks  we  shall  use  the  following  classical  expressions  for  the  momen¬ 
tum  and  energy  exchange  rates. 

For  the  momentum  exchange  rate  we  take  the  classical  Coulomb  (electron- 
ion)  collision  frequency[9]  for  momentum  relaxation  (i'?)^^ 


where  the  plasma  parameter  A  is  given  by 


(75) 


i  wnoAl 


This  collision  frequency  determines  the  classical  Spitzer  resistivity 


(76) 


noc* 


(77) 
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For  compatibility,  we  normalize  with  the  lower  hybrid  frequency  and  cast  the 
result  in  terms  of  our  dimensionless  parameters,  to  get 

^  _J_  ( miy^’^a>pelnA 
WM  \/2ir  Vm*/  Ua  A 

For  a  heating  rate  benchmark  we  define  a  classical  heating  rate. 

Joule  heating 

(‘'e  )ci  =  ^7;  3  noTe 

which  yields 

Wih  3\/^\VttJ  Tt\Tn^)  Uiet  A 

_  £  /  Ud«\  Tj  wig  (^'e)c/, 

3\vttJ  Term  Uih 

where  the  second  equation  shows  the  explicit  dependence  on  the  collision  fre¬ 
quency. 

Finally,  we  note  that  in  calculating  the  anomalous  rates  we  approximate 
the  integrals,  as  commonly  done  in  the  literature,  by  the  contribution  of  the 
dominant  mode  only  (t.e  for  fc**),  meaning  that  all  the  Fourier-decomposed 
properties  are  estimated  at  the  doubly  maximized  growth  (ie  maximized  over 
wavelength  and  propagation  angle)®. 

3.7.2  Numerical  Results 

Since  /3e  and  uje/vu  are  the  two  parameters  that  vary  most  within  the  plasma 
of  the  MPD  thruster,  they  were  chosen  as  the  varying  parameters  of  the  calcu¬ 
lations.  When  ^e  is  varied,  uje/vt,  is  kept  at  20,  and  when  uje/vi,  is  varied,  /3e  is 
set  at  unity.  The  other  parameters  are  rm/rrie  =  73300  (for  argon),  Ti/Te  =  1, 
t'e/w/fc  =  1  and  tjjpe/uee  =  100  for  continuity  with  the  calculations  in  of  our 
previous  studies’”. 

*AIthough  this  siinpIific4tion  is  followed  by  almost  all  anomalous  transport  calculations  in  the 
literature,  (often  without  justification)  one  must,  for  quantitative  accuracy,  check  to  see  whether, 
indeed,  the  Jb  resolved  rate  calculations  vary  weakly  with  the  wavenumber  in  the  region  of  parameter- 
space  of  interest.  Although  we  have  done  a  few  cursOTy  checks  on  the  sensitivity  of  the  rates  to  the 
details  of  the  spectrum  we  have  not  evaluated  the  above  calculation  method  over  a  comprehensive 
range  to  warrant  its  accuracy  over  the  entire  plasma  parameter  range  of  interest,  since  we  are  only 
concerned  here  with  order-of-magnitude  variations  in  the  rates.  More  accurate  calculations  in  the 
future  should  address  this  question. 

’Higher  collisional  rates  will  be  considered  in  the  MPD  thruster-specific  calculations  of  the  next 
section. 


(78) 

(79) 


(80) 
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Effects  of  Plasma  Beta  The  effects  of  plasma  beta  on  the  resistivity  are 
shown  in  Fig.  (2)  where  the  ratio  of  anomalous  to  classical  momentum  ex¬ 
change  frequency  (which  is  essentially  the  ratio  of  the  corresponding  resistivi¬ 
ties  (c/  Eqs.  (18)  and  Eq.  (77))  is  plotted  versus  beta  for  the  parameters  listed 
above.  We  immediately  note  that,  even  when  plotted  on  such  a  compacted  log- 


Figure  2:  The  anomalous  momentum  exchange  frequency,  normalized 

by  its  classical  counterpart  and  plotted  versus  the  electron  beta  according  to 
four  saturation  models.  Argon  with  =  20,  F./Ti  =  1,  Ve/^^ih  =  1  and 

^p€  / ^Ot  *  1 

arithmic  ordinate  scale  spanning  eight  orders  of  magnitude,  all  the  curves  can 
be  seen  to  deviate  from  their  *  0  asymptotes  (which  are  practically  reached 
at  =  .001).  These  deviations  are  due  to  the  electromagnetic  corrections  to 
the  electrostatic  limits,  as  separated  in  Eq.  (58).  We  also  note  that  the  trapping 
models  are  most  effected  by  finite-beta  effects. 

The  ion  trapping  model  is  of  special  interest  as  discussed  in  Section  3.6.1 
especially  since  it  was  the  only  one  assumed  in  the  purely  electrostatic  study 
of  ref.  [21].  We  see  that,  when  /?,  is  on  the  order  of  unity  or  greater,  as  is 
commonly  the  case  of  the  MPD  thruster  plasma,  the  anomalous  resistivity  is  an 
order  of  magnitude  less  than  that  predicted  by  the  purely  electrostatic  limit. 

The  reason  the  anomalous  resistivity  decreases  with  increasing  beta  accord¬ 
ing  to  trapping  models  can  be  traced  to  the  coupling  with  pseudo- whistler  modes 
that  we  discussed  in  ref.  [5].  As  beta  increases,  the  disturbances  to  the  mag¬ 
netic  field  do  not  have  the  time  to  dissipate  (low  Alfv^n  velocity)  and  significant 


CHOUEIRJ,  KELLY  and  JAHN 


22 


electromagnetic  coupling  arises.  The  unstable  waves  acquire  some  of  the  char¬ 
acteristics  of  oblique  whistlers  and  consequently  the  most  unstable  modes  shift 
to  lower  frequencies.  Since  the  saturation  level  due  to  trapping  scales  with  fre¬ 
quency  to  the  fourth  power  (c/  Eq.  (70)),  the  end  effect  is  a  substantial  reduction 
in  the  anomalous  resistivity. 

It  is  relevant,  in  the  context  of  a  comparative  study,  to  mention  thai  a 
plausibility  argument  based  on  the  concept  of  minimum  dissipation  is  usually 
invoked  in  the  literature  when  saturation  mechanisms  are  compared.  According 
to  this  argument  the  mechanism  requiring  the  least  energy  at  saturation  is 
favored.  A  comparative  analysis  based  solely  on  this  criterion  can,  however, 
be  misleading  if  it  does  not  consider  the  viability  of  the  physical  mechanism 
independent  of  the  global  thermodynamics. 

We  note  from  the  plot  that  the  Fowler  bound  on  the  calculated  rates  allows, 
in  principle,  for  a  wide  latitude  for  anomalous  resistivity  to  be  important. 

We  should  not  expect  the  electron  trapping  model  to  dictate  the  transport 
for  arguments  already  made  in  previous  sections.  Furthermore,  we  should  men¬ 
tion  that  more  careful  studies  of  resonance  broadening  than  was  made  at  the 
time  the  mechanism  was  first  proposed,  have  shown  that  its  effects  are  limited 
to  a  redistribution  of  energy  in  fc-space  at  low  plasma  beta,  and  that  it  does 
not  result  in  enough  dissipation  to  saturate  an  instability  (such  as  the  case  of 
the  lower  hybrid  gradient  driven  instability  (LHGDI),  for  instance)  .  Therefore, 
for  low  ion  trapping  is  the  most  viable  mechanism.  At  these  conditions, 
the  anomalous  resistivity  can  be  quite  dominant  (as  is  observed  in  Fig.  (2)), 
more  than  two  orders  of  magnitude  larger  than  the  classical  value,  in  agreement 
with  the  findings  of  ref.  [21].  As  beta  increases,  saturation  by  resonance  broad¬ 
ening  can  become  more  viable  especially  since  the  turbulent  saturation  levels 
are  considerably  lower  than  those  for  ion  trapping  (as  is  clear  from  the  same 
plot).  Whether  one  or  the  other  mechanism  controls  saturation  depends,  at  least 
partly,  on  whether  the  spectrum  is  narrow  or  broad.  Even  though  experimental 
data  on  turbulent  fluctuations  in  the  MPD  thruster[10,  11]  give  evidence  of  a 
dominant  narrow  (peaked)  spectrum  of  turbulence  in  the  lower  hybrid  range, 
the  considerably  lower  levels  of  saturation  energy  implied  by  the  ion  resonance 
broadening  mechanism  warrant  its  consideration  as  a  contender  in  the  control 
of  turbulent  transport.  Of  course,  this  question  is  best  answered  by  computer 
particle  simulations. 

If  ion  resonance  broadening  does  take  over  the  control  of  saturation,  anoma¬ 
lous  transport  can,  for  the  parameters  in  the  above  calculations,  be  brought 
down  below  classical  levels.  One  should  therefore  expect,  by  virtue  of  the  sub¬ 
stantial  variability  of  the  plasma  beta  within  the  MPD  thruster  plasma,  that 
there  are  regions  where  anomalous  resistivity  dominates  over  its  classics'  analog 
as  there  are  regions  where  the  converse  is  true.  More  on  this  issue  will  be  said 
in  Section  3.7.3. 

The  same  comments  we  made  above  also  apply  for  the  anomalous  (  ec- 
tron  heating  rates  that  were  normalized  by  the  joule  heating  rate  and  plotted 
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in  Fig.  (3).  It  is  clear  from  this  plot  that  when  ion  trapping  dominates,  the 


Figure  3:  The  anomalous  electron  heating  rate,  normalized  by  the  joule 

healing  rate  and  plotted  versus  the  electron  beta  according  to  four  saturation 
models.  Argon  with  =  20,  T./T,  =  1,  =  1  and  ufpefufce  =  100. 

anomalous  heating  rate  is  substantially  larger  than  the  electron  joule  heating 
rate.  Ion  heating  rates  are  not  shown  here  but  were  calculated  in  ref.  [5]  and 
were  found  to  be  similar  in  both  magnitude  and  dependence  as  their  electron 
counterparts. 

To  compare  the  two  rates  we  have  calculated  their  ratio  and  plotted  the 
result  in  Fig.  (4).  There  is  only  one  curve  in  this  figure  because  the  various 
saturation  models  cancel  out  in  the  division.  Since  this  ratio  is  independent 
of  the  saturation  details,  it  is  more  accurate  than  the  other  quantities  we  have 
calculated.  We  note  from  this  figure  that,  in  the  electrostatic  limit,  the  two 
anomalous  heating  rates  are  basically  equal.  This  feature  is  in  contrast  to  the 
way  electrons  and  ions  are  heated  classically  (especially  for  a  heavy  atom  like 
argon)  and  is  a  well-known  characteristic  of  the  electrostatic  LHCDI  (or  the 
modified  two-stream  instability,  MTSI)  as  noted  in  ref.  [22,  23).  Since  the  ions, 
are  heated  by  the  instability-induced  turbulence  at  rates  comparable  to  those 
of  the  electrons,  and  since  in  the  MPD  thruster,  the  electron  energy  is  strongly 
tied  to  excitation  and  ionization  through  inelastic  collisions,  anomalous  heating 
can  offer  an  answer  to  the  long  standing  question  of  why  the  ion  temperature  is 
somewhat  higher  than  the  electron  temperature.  Of  course,  for  this  argument 
to  be  true  not  only  {I'D  must  be  comparable  to  {I'D  an  '  saturation 
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Figure  4:  The  anomalous  electron  heating  rate,  normalized  by  the 

anomalous  ion  heating  rate,  plotted  versus  the  electron  beta.  Argon 

with  Uie/vti  =  20,  Ti/Tt  =  1,  Ve/uiih  =  1  and  =  100. 


level  must  be  high  enough  to  warrant  the  dominance  of  anomalous  heating  over 
classical  heating.  Such  is  the  case  when  the  instability  saturates  by  trapping 
ions. 

The  above  argument  about  the  relative  temperatures  is  strongest  in  the 
electrostatic  limit  and  is  in  agreement  with  ref.  [21].  When  electromagnetic 
effects  start  to  become  important  with  increasing  beta,  the  same  figure  shows 
a  degradation  of  the  heating  parity  towards  a  progressively  preferential  heating 
of  electrons.  This  finding  is  in  agreement  with  that  of  ref.  [16]  where  only  the 
collisionless  limit  was  studied.  This  degradation  in  heating  parity  is  not  strong 
enough,  however,  to  weaken  the  grounds  for  the  above  argument  concerning  the 
relative  temperatures,  especially  for  a  heavy  atom  like  argon.  Indeed,  we  see 
from  the  same  figure  that  a  four-octave  increase  in  beta  does  not  change  the 
order  of  magnitude  of  the  relative  heating  ratio. 

The  increase  of  preferential  electron  heating  with  increasing  beta  may  be 
partly  due  to  the  fact  that,  at  low  beta,  the  instability  has  its  dominant  modes 
oriented  at  small  angles  to  the  magnetic  field  (kg/k  ^  or  cs  1)  as 

was  found  in  our  previous  studies[3],  and  consequently  “perceives”  the  electron 
with  an  effective  mass  comparable  to  that  of  the  ions[22].  As  beta  increases,  elec¬ 
tromagnetic  coupling  with  oblique  pseudo-whistlers  cause  the  dominant  modes 
to  propagate  more  obliquely,  as  first  noted  by  refs.  [24]  and  [25]  and  extended 
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by  us,  in  ref.  [5],  to  the  collisional  domain.  Consequently,  the  effective  electron 
mass  decreases  and  the  electrons  become  much  easier  to  heat  than  the  ions. 

Effects  of  the  Drift  Velocity  The  effects  of  the  drift  velocity  are  illustrated 
in  the  plots  of  Figs.  (5)  and  (6)  for  the  same  parameters  as  above  but  with  Pe 
set  to  unity  and  Ude/ut,  varying  between  10  and  100. 


“de/Vti 


Figure  5:  The  anomalous  momentum  exchange  frequency,  normalized 

by  its  classical  counterpart  and  plotted  versus  the  normalized  drift  velocity 
according  to  four  saturation  models.  Argon  with  Pe  =  1.  Ti/Te  =  1,  Vef^ih  =  1 

and  <x>pel<^ae  =  100. 

In  reference  to  Fig.  (5)  we  note  that  the  general  decreasing  trend  of  anoma¬ 
lous  resistivity  with  the  drift  velocity,  once  the  instability  is  onset®,  is  not  in¬ 
tuitive.  One  would  expect  that  an  increase  in  the  free  energy  source  of  the 
instability  would  enhance  the  anomalous  resistivity  effect.  In  ref.  [21]  the  same 
trend  was  found  but  no  explanation  was  given. 

This  trend  can  be  understood  once  we  realize  that  the  scaling  of  the  linear 
growth  rate  of  the  dominant  mode  (whidi  does  increase  with  the  drift  velocity) 
does  not  necessarily  reflect  in  weak  turbulence  (quasi-linear)  transport  scaling 
since  the  dependencies  of  the  saturation  mechanism  (which  is  extraneous  to  lin¬ 
ear  theory)  can  overwhelm  linear  trends.  This  becomes  clearer  by  noting  that 
although  an  increase  in  the  drift  velocity  does  enhance  the  linear  growth  rate 

*The  u^/vti  thresholds  for  the  onset  of  the  insubility  (start  of  positive  linear  growth)  are  not 
marked  on  these  plots  because  they  are  on  the  order  of  unity. 
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of  the  dominant  mode  as  found  in  a  previous  study  of  ours{3],  it  also  shifts  the 
modes  to  more  oblique  propagation  and  lowers  their  frequencies.  Even  though 
the  instability  goes  to  longer  wavelength[5],  the  dependence  of  the  saturation 
level  for  a  trapping  mechanism  (c/  Eq.  (70))  scales  with  the  frequency  to  the 
fourth  power  so  that  the  frequency  scaling  of  the  saturation  mechanism  over¬ 
powers  the  growth  scaling  of  the  linear  modes®. 

The  Fowler  bound,  on  the  other  hand,  whose  scaling  is  more  strongly  tied 
to  the  free  energy  source,  does  increase  monotonically  with  the  drift  velocity,  as 
seen  from  the  resistivity  plot. 

As  expected  from  the  above  beta-dependence  study,  anomalous  electron 
heating  rates  exceed  those  of  the  ions  for  the  present  case  of  =  1-  The 
preferential  electron  heating  is  further  enhanced  by  increasing  drift  velocity  as 
can  be  seen  in  Fig.  (6).  The  reason  for  this  behavior  is  similar  to  the  one  given 


Figure  6:  The  anomalous  electron  heating  rate,  normalized  by  the 

anomalous  ion  heating  rate,  and  plotted  for  the  same  conditions  as  in  Fig.  (5). 

above  in  the  context  of  electromagnetic  enhancement  of  electron  preferential 
heating.  This  is  so  because  both  increasing  beta  and  increasing  drift  velocity 
act  to  shift  the  instability  toward  more  oblique  propagation  thus  reducing  the 

*This  trend  is  further  accentuated  for  saturation  by  dectron  trapping  because,  in  addition  to  the 
above  arguments,  the  saturation  level  scales  with  and  4’  decreases  considerably  (oblique  propa¬ 
gation)  with  increasing  drift  velocity.  At  very  high  drift,  the  Doppler  shift  term  in  the  saturation 
model  (c/  Eq.  (72))  becomes  more  significant  and  reverts  the  trend,  which  explains  the  rise  of  the 
electron  trapping  curve  in  Fig.  (5)  at  high  values  of  usi/sfa. 
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large  if»e-effect  (where  m*  is  the  effective  electron  mass  that  scales  with  the 
square  of  V’)  and  subjecting  the  now  “lighter”  electrons  to  more  heating. 

3.7.3  MPD  Thruster  Calculations 

We  have,  in  the  above  calculations,  chosen  a  set  of  parameters  that  is  generally 
representative  of  the  MPD  thruster  plasma.  There  is,  however,  one  exception.  It 
is  the  value  of  Vt/uih  which  we  have  set  equal  to  unity  as  a  compromise  between 
having  to  represent  a  collisional  plasma  and  providing  a  link  with  previous  stud¬ 
ies.  Moreover,  many  of  the  complex  interactions  between  the  natural  plasma 
modes,  the  free  energy  source  and  collisions  are  most  pronounced  when  the  colli¬ 
sion  frequency  is  on  the  order  of  the  oscillating  frequency.  We  now,  supplement 
our  calculations  with  results  obtained  at  collisional  levels  more  appropriate  of 
the  MPD  thruster  plasma. 

In  order  to  approximate  a  typical  range  for  MPD  thruster  plasma  colli- 
sionality,  we  consider  the  typical  range  for  the  variation  of  temperature  and 
density.  For  more  detail  on  how  the  various  parameters  of  interest  vary  within 
the  MPD  thruster  discharge  the  reader  is  referred  to  the  recent  parameter  re¬ 
view  in  ref.  [10].  Assuming  that  T*  varies  between  1.5  and  3  eV,  while  no  ranges 
between  10“  and  1.5  x  10^  m"^,  we  can  calculate  a  lower  and  upper  bound 
for  i^e/^tk  in  argon  from  Eq.  (78)  to  be  25  and  500,  respectively,  where  we  have 
fixed  LjftfuJa!  at  100  for  compatibility  with  the  above  calculations. 

For  the  results  shown  in  Fig.  (7)  we  have  chosen  to  fix  beta  at  unity  to 
preserve  electromagnetic  effects  and  varied  uoe/vti  from  100  down  to  the  thresh¬ 
old  of  the  instability,  which,  although  slightly  exaggerated  in  the  figure,  was  at 
ti* /va  2;  1 .5.  For  each  of  the  three  considered  mechanisms  the  plot  shows  a  band 
whose  upper  line  corresponds  to  the  moderately  collisional  condHion  =  25 

(Tg  =  3  eV,  no  =  10^  m"^)  and  whose  lower  (broken)  line  corresponds  to  the 
strongly  collisional  condition  v^fujih  =  500  (T*  =  1.5  eV,  no  =  1.5  x  10®  m”^). 

The  trends  in  this  plot  are  similar  to  those  already  illustrated  in  Fig.  (5) 
which  indicated  that  we  have  not  jeopardized  any  physics  by  assuming  v^luiih  — 
1  previously. 

We  note  from  the  figure  that,  although  the  Fowler  bound  allows  for  a  large 
microturbulent  contribution  to  the  resistivity,  ion  resonance  broadening  might 
cause  the  instability  to  saturate  at  low  levels.  Even  though  arguments  have  been 
advanced  recently  discounting  the  effectivity  of  sudi  a  mechanism,  it  should  not 
be  totally  discounted  pending  strong  evidence  from  computer  particle  simula¬ 
tions  and/or  dedicated  experiments. 

We  furthermore  see  that,  in  the  case  of  ion  trapping  saturation,  once  the 
instability  is  onset,  the  importance  of  anomalous  resistivity  in  the  MPD  thruster 
plasma  is  not  as  much  dictated  by  the  drift  velocity  (since  the  two  limiting  curves 
are  quite  fiat),  as  one  would  intuitively  suspect,  as  it  is  dictated  by  the  level  of 
collisionality.  Indeed,  if  collisionality  is  strong,  anomalous  resistivity  can  be  kept 
below  classical  levels  even  if  the  instability  is  excited  and  even  if  ion  trapping 
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Figure  7:  The  anomalous  momentum  exchange  frequency,  normalized  by 

its  classical  counterpart  and  plotted  versus  the  normalized  electron  drift  velocity 
according  to  three  saturation  models.  The  upper  line  of  eadi  band  corresponds 
to  the  moderately  collisional  condition  =  25  (T*  =  3  eV,  no  =  10“  m"^) 
and  the  lower  (broken)  line  corresponds  to  the  strongly  collisional  condition 
z/./wjh  =  500  (Te  =  1.5  eV,  no  =  1.5  X  lO*^  m'^).  Argon  with  ^,  =  1,  Ti/Te  =  I, 
and  K  100. 


is  responsible  for  saturation,  as  is  clear  from  the  plot“.  This  implies  that  low 
density  regions  of  the  MPD  thruster  discharge,  such  as  regions  depleted  from 
charge  due  to  the  Lorentz  force  component,  tend  to  be  more  vulnerable 
to  anomalous  resistivity  than  denser  (or  more  collisionally  dominated)  regions. 
This  trend  is  in  agreement  with  the  well-known  fact  that  dissipation  in  charged- 
dcpleted  regions  of  the  device,  like  the  anode  vicinity[26],  is  enhanced  by  weak 
collisionality. 

Stated  differently,  under  MPD  thruster  plasma  conditions  and  for  the  mi¬ 
croinstabilities  in  question,  the  level  of  anomalous  contribution  to  resistivity  is 
dictated  to  a  large  extent  by  the  parameter  It  is  interesting  to  note  that 

this  parameter  is  directly  related  to  the  electron  Hall  parameter  D//,.  Indeed, 
it  is  just  the  inverse  of  the  electron  Hall  parameter  scaled  by  the  square  root  of 
the  mass  ratio 


O  _  Woe  _ 

~  1/.  ~  i/./wu 


(81) 


"’ll  must  be  said,  however,  that  even  in  the  case  of  high  coUisionality  ~  500)  where 

anomalous  resistivity  is  kept  below  classical  levels,  it  is  still  a  6nite  fraction  of  its  classical  counter¬ 
part  (about  25%  in  the  above  calculations),  as  can  be  seen  from  the  same  plot  (again,  assuming  ion 
trapping  saturation). 
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The  known  scaling  of  the  anode  voltage  drop  with  the  Hall  parameter  (see 
for  instance  the  recent  measurements  in  ref.  (26])  that  constitutes  one  of  the 
most  dissipative  sinks  for  the  low-power  MPD  thruster  is  thus  another  invariant 
behavioral  trait  of  the  accelerator  that  could  possibly  be  explained  by  the  effects 
of  microinstabilities. 

The  anomalous  ion  and  electron  heating  rates  were  also  calculated  in  ref.  [5] 
and  found  to  have  the  same  general  trends  as  those  of  the  anomalous  resistivity. 

We  now  have  all  that  we  need  to  cast  the  anomalous  transport  models  in 
a  form  that  can  be  used  in  plasma  fluid  codes.  This  will  be  documented  in  th 
enext  sections. 


4  Anomalous  Transport  Models  and  their  In¬ 
clusion  in  Fluid  Codes 

4.1  Goals 

A  spatial  resolution  of  the  realms  of  anomalous  and  classical  transport  within 
the  MPD  thruster  discharge  and  their  dependence  on  the  accelerator’s  operating 
conditions  can  only  be  addressed  accurately  with  an  improved  fluid  code  that  in¬ 
corporates  the  above  theories  in  a  self-consistent  fluid-kinetic  description.  This 
may  be  accomplished  by  carrying  a  priori  calculations  of  the  relevant  anoma¬ 
lous  transport  for  the  expected  parameter-space  covered  by  typical  numerical 
simulations  then  fitting  the  calculations  with  polynomial  expressions.  These 
expressions  become  the  transport  models  suitable  for  inclusion  in  codes  for  the 
self-consistent  numerical  simulation  of  MPD  thruster  flows. 


4.2  Final  Models 


In  general  the  microstability  (and  hence  microturbulence)  description,  detailed 
in  the  sections  above,  depends  on  the  following  set  of  eight  independent  macro¬ 
scopic  parameters 


Uh. 

m«  ’ 


> 


(82) 


The  first  two  parameters  kr„  (r„  being  the  electron  cyclotron  radius)  and  4^ 
represent  the  normalized  wavenumber  and  propagation  angle  (with  respect  to 
the  magnetic  field)  of  the  oscillations  respectively  and  are  varied  to  growth- 
maximize  the  solutions.  Since  all  anomalous  transport  rates  used  here  were 
calculated  at  maximum  growth  these  two  parameters  drop  out  of  the  final  mod¬ 
els.  The  mass  ratio  nii/me  is  that  of  argon.  AH  solutions  were  found  in  the 
calculations  of  the  previous  sections  to  be  very  insensitive  to  the  fourth  param¬ 
eter,  namely  the  ratio  of  electron  plasma  frequency  to  the  electron  cyclotron 
frequency  Wpt/uo,,  as  long  as  that  ratio  exceeded  10  which  was  the  case  for  the 
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simulations  conducted  so  far[27,  28].  Similarly,  the  solutions  were  weakly  depen¬ 
dent  on  Pc  (the  ratio  of  electron  thermal  pressure  to  magnetic  pressure)  as  long 
as  the  electron  Hall  parameter  did  not  exceed  10.  Although  that  was  the  case 
for  the  numerical  simulations  in  Refs.  [27,  28],  it  is  expected  that  the  simulation 
of  more  realistic  geometries  at  high  total  currents  would  raise  the  electron  Hall 
parameter  enough  to  require  the  full  inclusion  of  finite-beta  effects. 

The  last  three  parameters  are  the  most  important  for  our  problem.  First, 
uje/vti  must  reach  a  threshold  for  the  instability  to  be  excited  and  hence  for 
anomalous  transport  to  be  operative.  For  the  entire  region  of  the  investigated 
parameter-space  that  threshold  was  very  near  1.5.  Second,  the  ion  to  electron 
temperature  ratio  plays  a  role  in  scaling  the  level  of  turbulence.  Invariably  for 
our  parameter-space,  it  was  found  that  increasing  Ti/Tc  causes  a  devaluation  of 
anomalous  transport.  The  most  important  of  all  the  macroscopic  parameters 
turned  out  to  be  the  last  one  namely  the  electron  Hall  parameter  Hh.- 

The  anomalous  resistivity  t)an 


Van  = 


^c{nc)an 
(?nc  ' 


(83) 


calculated  using  the  theory  presented  in  the  above  sections,  and  normalized  by 
its  classic  counterpart  vci  =  Tncud^n,  is  shown  in  Fig.  (8).  It  is  important  to 
note  that  an  increase  in  the  electron  Hall  parameter  for  typical  values  of  Ti/Tt 
leads  to  a  very  significant  increase  in  the  anomalous  resistivity  if  the  parameter 
tide/vu  is  above  the  stability  threshold.  It  is  interesting  to  note  that  the  scaling 
of  this  ratio  with  the  Hall  parameter  is  in  general  agreement  with  that  recently 
inferred  by  Gallimore|29]  from  measurements  in  the  anode  region. 

A  similar  plot  is  shown  in  Fig.  (9)  for  the  ion  heating  rate  {uI)an  normalized 
by  the  Coulomb  frequency. 

A  two-parameter,  variable  cross-term,  least  square  fit  was  made  to  the  cal¬ 
culated  rates  shown  in  Figs.  (8)  and  (9)  in  order  to  make  them  suitable  for 
inclusion  in  plasma  fluid  flow  codes. 

The  resulting  two-parameter  interpolating  polynomial  for  («'7^)xw/t'ei  has  an 
average  accuracy  of  15%  and  reads 


=  5.36  X  10-®  +  1.29  X  10-*  Q 

^ei 

+  6.03  X  10-«  fl**  +  9.44  X  lO”*  fi® 

+  ^(-7.55  X  lO-'^  -  5.41  X  10-*  fl 

-  3.93  X  10-*  n*).  (84) 

The  ions  are  heated  by  the  turbulent  fluctuations  at  a  rate  (Qi)x^  =  §(»'7^)>«wT>, 
The  effective  conductivity  introducing  the  anomalous  resistivity  effect  to  the 
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Figure  8:  Ratio  of  anomalous  resistivity  to  classical  resistivity  as  a  function  of 
the  electron  Hall  parameter  and  Ti/T*  with  uje/va  exceeding  1.5, 


flow  code  has  the  form 


V/  = 


C^Tle 


me{l'ei  +  {u^)ASy 


(85) 


where  is  the  electron-wave  momentum  exchange  frequency,  which  is 

again  computed  through  an  interpolating  polynomial  of  average  accuracy  of 


10% 


I'ei 


0.192  -f-  3.33  X  10-2  Q  ^  21202 

-  8.27  X  10"*  n®  +  ^(1.23  X  10-2 

■*e 

-  1.58  X  10-2  Q 

-  7.89  X  10-2  q2j 


(86) 


The  use  of  these  models  in  fluid  code  may  proceed  in  the  following  way.  At 
all  the  gridpoints  where  «*/««  <1.5  both,  \vg)AN  and  {u[)as  are  set  to  zero 
and  all  transport  is  assumed  purely  classical.  Otherwise,  the  anomalous  rates 
are  computed  from  the  above  polynomials  using  the  instantaneous  macroscopic 
parameters  and  folded  back  into  the  flow  equations  at  every  time  step  thus 
insuring  self-consistency. 
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Electron  Hall  Parameter 


Figure  9:  Anomalous  ion  heating  rate  normalized  by  the  Coulomb  frequency  as 
a  function  of  the  electron  Hall  and  7i/7i  with  ujc/vu  exceeding  1 .5 

4.3  Fluid  Code  Description 

The  plasma  fluid  code  used  to  test  the  anomalous  transport  models  derived 
above  was  developed  at  our  laboratory  partly  under  the  AFOSR  grant.  It  is 
described  in  more  detail  in  refs.  [27,  28].  We  provide  here  a  brief  description 
of  the  code.  Our  focus  will  instead  be  on  the  calculations  made  with  that  code 
under  the  support  of  AFOSR. 

The  results  presented  in  this  final  report  are  for  a  cylindrical  MPD  thruster- 
geometry.  The  present  code,  however,  can  handle  any  axi-symmetric  geometry. 

The  flow  field  code  uses  a  finite-volume  discretization  with  artificial  dissi¬ 
pation  described  by  Jameson|30].  The  time  stepping  is  done  via  the  explicit 
Euler  Forward  method,  and  convergence  is  accelerated  using  a  multigrid  itera¬ 
tion  first  proposed  by  Jameson  and  Jayram|31].  These  methods  yield  a  second 
order  steady-state  solution  to  the  conservation  equations. 

The  solution  of  the  conservation  equations  through  the  finite  volumes/Euler 
forward  routine  is  alternated  with  the  solution  of  the  electromagnetic  equation, 
until  consistency  among  all  the  parameters  is  achieved. 

The  stream  function  equation  is  solved  with  a  second-order  nonlinear  ex¬ 
plicit  scheme  develop>ed  at  EPPDyL.  After  a  central  difference  discretization  is 
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performed,  the  stream  function  equation  assumes  the  form 

Cl^y  +  <2  =  0.  (87) 

Here  Vy  stands  for  the  value  of  ip  at  the  horizontal  and  j**  vertical  discrete 
steps.  The  coefficients  Ci  and  then,  are  functions  of  the  differential  equation’s 
coefficients  and  of  V’i.j+i  Ea<^  iteration  step  we  set 

(88) 

until  (V'y)^"'^'V(V'y)^”^  is  equal  to  1.  Using  as  initial  guess  the  value  of  ip  com¬ 
puted  at  the  previous  MacCormack  step,  convergence  is  achieved  in  a  number 
of  steps  ranging  between  1  and  250. 

These  methods  represent  an  improvement  with  respect  to  the  MacCormack 
method  used  previously  since  it  decreases  running  time  by  up  to  a  factor  of  eight 
while  reducing  the  steady-stale  residual.  The  model  was  coded  in  APL2  (taking 
full  advantage  of  the  inherent  vectorizing  capability  of  the  language)  and  has 
been  run  on  a  variety  of  machines  ranging  fi-om  a  Mac  Quadra  (1.2  Mflops)  to 
an  IBM  ES/3090  600J  Supercomputer  where  specialized  APL2  compilers  can 
allow  convergence  to  be  reached  in  less  than  30  minutes  of  CPU  time. 

4.4  Results 

This  section  will  first  discuss  the  performance  curves  obtained  running  at  various 
current  levels,  with  and  without  anomalous  transport,  and  compare  them  to 
experimental  data  taken  by  Burton  for  a  FSBT(32].  Two  high  current  runs 
(18  kA),  one  with  classical  and  one  with  anomalous  dissipation,  will  then  be 
compared,  in  order  to  asses  in  what  areas  of  the  thruster  anomalous  transport 
has  the  greatest  impact.  Furthermore,  current  distribution  patterns  will  be 
compared  with  measured  ones.  The  current  density  distribution  for  simulations 
with  and  without  anomalous  transport  will  also  be  shown,  together  with  an 
experimental  distribution  obtained  by  Cilland|33]. 

Henceforth,  the  runs  in  which  the  conservation  equations  include  anoma¬ 
lous  transport  coefficients  will  be  denoted  as  “anomalous”,  whereas  those  not 
containing  anomalous  transport  will  be  called  “classical  runs" . 

4.4.1  Discussion  of  Performance  Curves 

The  code  was  run  with  a  mass  flow  rate  of  6  g/s  and  with  currents  ranging 
from  7  kA  to  18  kA.  The  code  diverged  for  higher  currents.  Runs  were  made 
both  with  and  without  anomalous  transport  dissipation  terms  present  in  the 
calculation.  Plots  of  thrust,  voltage  and  efficiency  versus  J^/rh  are  shown  for 
these  runs  along  with  experimental  data  when  appropriate. 


I 


Figure  10;  Calculated  thrust  from  simulations  with  and  without  anomalous 
transport  compared  with  measured  FSBT  thrust  measured  by  Gilland[33]  and 
Burton  [32]. 


As  is  evident  from  Fig.  (10),  thrust  scales  as  the  square  of  the  current  in¬ 
dividually  for  the  case  with  and  without  anomalous  transport.  Both  classical 
and  anomalous  predictions  compare  very  well  with  the  measured  Benchmark 
values,  which  in  turn  are  very  close  to  the  theoretical  Maecker  thrust|34]  for 
the  higher  range  of  specific  impulse  (J^/m  >  40  kA®  sec/g).  At  lower  values 
of  J^/rh,  the  measured  thrust  data  exceed  the  Maecker  law  prediction.  The 
simulations,  h^'wever,  were  able  to  account  for  that  excess  as  can  be  seen  from 
the  figure.  It  is  also  important  to  note  that  the  thrust  from  the  anomalous  run 
becomes  appreciably  lower  as  J^/rh  increases.  This  may  be  attributed  to  the 
fact  that  for  higher  Tj,  (anomalous  case),  the  gas  pressure  decreases  less  along 
the  chamber,  so  that  the  -Vp*  force  term  is  reduced. 

Figure  (11)  shows  the  calculated  voltage  drop  across  the  plasma  (from  the 
anode  to  the  cathode).  The  calculated  voltages  do  not  include  electrode  sheath 
falls  and  therefore  cannot  be  easily  compared  to  the  measured  terminal  voltage 
of  the  thruster.  The  voltage  starts  being  consistently  higher  in  the  anomalous 
case  only  for  the  higher  currents  (above  around  40  kA*  sec/g). 

The  increase  in  the  voltage  and  the  decrease  in  thrust,  due  to  anomalous 
transport,  translate  directly  into  a  degradation  in  the  thrust  efficiency.  To  see 
the  extent  of  this  degradation,  the  thrust  efficiency  is  calculated  through  the 


formula 


F* 

2m JV 


(89) 


In  Fig.  (12)  the  efficiency  is  plotted,  again,  versus  J®/m  .  While  for  currents 
below  10  kA  (J*/rh  <  20  kA*  sec/g)  there  is  hardly  any  difference  between  the 
classical  and  the  anomalous  cases,  for  higher  currents  the  efficiency  is  substan¬ 
tially  lower  with  anomalous  transport  dropping  by  13%  at  J*/m  =  55  kA*  sec/g. 
Furthermore,  efficiency  from  anomalous  runs  tends  to  reach  a  plateau  around 
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Figure  11:  Plasma  fall  for  simulations  with  and  without  anomalous  transport. 
The  calculated  voltage  does  not  include  electrode  sheath  falls. 


Figure  12:  Thrust  efficiency  for  simulations  with  and  without  anomalous  trans¬ 
port  compared  with  measured  FSBT  values. 
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Figure  13;  Current  contour  lines,  in  percentage  of  enclosed  total  current,  as 
measured  by  Gilland[33].  (J*/rh  =  30  kA^  sec/g). 

25%  for  J*/ih  >  30  kA^  sec/g,  while  in  the  classical  transport  case  it  keeps 
increasing.  The  shape  of  the  efficiency  curve  for  the  anomalous  runs,  then, 
is  closer  to  the  plateau  shape  of  the  experimental  curve.  Even  the  lower  effi¬ 
ciencies  predicted  by  the  simulation  with  anomalous  transport  are  higher  than 
the  experimentally  measured  values.  This  discrepancy  is  due  to  a  large  extent 
to  the  fact  that  this  numerical  simulation  does  not  include  the  voltage  drops 
from  the  electrode  sheaths.  If  the  experimental  total  voltage  were  reduced 
by  the  sheath  voltage,  the  numerical  and  experimental  curves  would  approach 
each  other  significantly.  Gallimore[29],  in  fact,  has  shown  that  the  anode  fall 
voltage  increases  monotonically  with  J*/rh  ,  and  reaches  values  as  high  as  40 
volts  around  J^/rh  =  50  kA*  sec/g.  Greater  accuracy  may  be  expected  with 
the  use  of  a  more  accurate  geometric  thnister  configuration  (with  a  varying 
cross-section)  through  the  introduction  of  a  transformation  of  coordinates  in 
the  model’s  equations.  This  will  be  the  next  step  in  our  future  studies.  Further¬ 
more,  the  extent  to  which  viscous  effects  may  help  approaching  the  experimental 
values  will  also  be  explored  in  future  investigations. 

4.4.2  Current  Distribution 

Current  distribution  in  a  full-scale  Princeton  Benchmark  Thruster  for  J®/ih  = 
30  kA®  sec/g  is  shown  in  Fig.  (13)  as  measured  by  Gilland[33j.  Current  dis¬ 
tribution  as  predicted  by  the  anomalous  transport  numerical  simulation  is  also 
shown  (Fig.  (14)). 

It  can  be  noted  from  these  two  figures  that,  along  the  cathode,  which  is  the 
region  where  the  geometric  configurations  are  most  similar,  the  calculated  cur¬ 
rent  distribution  is  evenly  spread  as  in  the  experiments,  with  attachment  going 
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Figure  14:  Current  contour  lines,  in  percentage  of  enclosed  total  current,  as 
predicted  by  a  flow  model  with  anomalous  transport.  (J*/m  =  30  kA*  sec/g). 


all  the  way  to  the  root.  Furthermore,  the  current  pattern  inaide  the  discharge 
chamber  is  well  predicted  both  in  shape  (contour  lines  curvature  and  inclina¬ 
tion)  and  magnitude.  This  is  the  case  up  to  the  60%  line  downstream  of  which 
the  calculated  current  does  not  “balloon"  out  as  much  as  in  the  experiments 
undoubtedly  because  of  the  presence  of  the  downstream  insulators  imposed  on 
the  modeled  geometry.  Particularly,  the  insulator  downstream  of  the  anode  does 
not  allow  attachment  to  the  anode  front  like  in  the  experiments. 

4.4.3  Comparison  Between  a  Classical  and  an  Anomalous  Run  for 
High  Current 

We  will  now  present  two  runs,  one  with  classical  and  one  with  anomalous  trans¬ 
port,  at  a  current  level  of  18  kA  and  a  mass  flow  rate  of  6  g/s,  and  discuss  the 
most  striking  differences  between  the  two.  In  particular,  resistivity  and  plasma 
heating  will  be  shown  to  be  significantly  enhanced  by  plasma  microturbulence. 

Of  course,  changing  the  conductivity  as  a  result  of  anomalous  effects  may 
affect  the  current  attachment  on  the  electrodes.  Figures  (15)  and  (16)  show  the 
distribution  of  current  for  the  classical  and  for  the  anomalous  runs.  Since  the 
conductivity  is  higher  in  the  plume,  for  the  classical  case  the  current  is  blown 
more  downstream.  Furthermore,  there  is  greater  attadiment  at  the  tip  of  the 
anode  for  the  classical  run,  again  because  of  the  higher  relative  conductivity 
there  as  will  become  apparent  from  the  resistivity  plot  shown  further  below. 

As  mentioned  above,  anomalous  transport  is  possible  only  in  regions  where 
Uit/vti  exceeds  1.5.  Moreover,  once  anomalous  transport  effects  are  onset  in  a 
certain  region,  their  magnitude  increases  monotonicaJly  with  the  electron  Hall 
parameter  as  can  be  seen  from  Eqs.  (8)  and  (9).  Figure  (17)  shows  a  map  of  the 
electron  Hall  parameter,  (!«,  in  locations  at  which  anomalous  transport  arises 
(i.e.  regions  with  Ude/vtt  >  1-5).  Fbr  such  high  a  current  as  18  kA,  and  unlike 
the  case  investigated  in  our  last  study[27],  anomalous  transport  can  be  expected 
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Figure  15:  Enclosed  current  lines  for  the  simulation  with  classical  transport 
(J  =  18  kA,  m  =  6  g/s).  The  number  along  each  line  indicates  what  percentage 
of  the  total  current  is  upstream  of  the  line. 


Figure  16:  Enclosed  current  lines  for  the  simulation  with  anomalous  transport 
(J  =  18  kA,  m  =  6  g/s). 
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to  affect  a  substantial  portion  of  the  discharge.  The  plot  in  Fig.  (17)  directly 
allows  the  identification  of  critical  regions  of  the  discharge  from  the  point  of 
view  of  anomalous  transport.  It  is  evident  that  the  cathode  root  and  tip  and 
the  anode  tip  are  most  critical. 

The  criticality  of  these  regions  is  further  supported  by  the  resistivity  map 
of  Fig.  (18)  where  the  ratio  of  “anomalous”  to  “classical”  resistivities  (i.e.  us¬ 
ing  the  resistivities  calculated  from  the  anomalous  and  classical  runs)  is  shown. 
The  pronounced  anomalous  resistivity  enhancement  in  the  immediate  vicinity 
of  the  anode  tip  is  very  reminiscent  of  the  enhanced  resistivities  inferred  by 
Gallimore[29]  from  local  measurements  near  the  anode.  This  further  substanti¬ 
ates  our  earlier  speculations[5, 27]  that  a  substantial  part  of  the  infamous  “anode 
drop”  may  be  related  to  the  turbulent  effects  induced  by  microinstabilities. 

An  appredation  of  the  extent  and  regional  preference  of  anomalous  plasma 
heating  can  be  obtained  by  comparing  the  maps  of  Fig.  (19)  and  Fig.  (20).  The 
two  maps  are  of  the  ion  temperature  from  the  classical  and  anomalous  runs 
respectively.  It  is  clear  from  these  figures  that  anomalous  transport  can  cause 
enhanced  plasma  heating  in  spedfic  regions.  The  “hot  spot”  near  the  badcplate 
at  the  cathode  root  has  been  observed  in  our  previous  study[27]  and  has  been 
suspected  to  play  a  role  in  backplate  erosion.  Badcplate  erosion  is  the  only 
erosion  mechanism  that  has  been  shown  experimentally  [35,  36]  to  occur  along 
with  the  appearance  of  “instabilities”  In  the  terminal  voltage  of  the  thruster. 
It  can  be  noted  from  Fig.  (19)  that  in  the  absence  of  microinstabilities  this 
“hot-spot”  disappears.  The  presence  of  such  a  spot  has  not  been  established 
experimentally  yet  and  could  provide  a  check  on  the  predictive  ability  of  the 
model. 

Finally,  we  show  in  Fig.  (21)  (7</7i  from  the  classical  run)  and  Fig.  (22) 
{Ti/Te  from  the  anomalous  run)  that  microturbulence  may  be  at  least  partly 
responsible  for  the  fact  that  7i/7i  can  exceed  unity  in  the  MPD  thruster  as  is 
well  known  from  experiments.  This  ratio  never  exceeded  unity  in  our  classical 
runs. 

We  have  demonstrated  that  performance  predictions  from  models  with  mi- 
croturbulent  transport  compare  more  favorably  with  experiments  than  do  pre¬ 
dictions  from  models  without  such  anomalous  effects.  Furthermore,  the  distri¬ 
bution  of  current,  resistivity  and  heavy  species  temperature  within  the  thrust 
chamber  shows  heat  dissipation  patterns  which  may  explain  some  features  of 
MPD  thruster  of>eration  such  as  the  anode  drop  and  backplate  erosion. 

Among  the  improvements  that  can  be  applied  to  the  numerical  model  pre¬ 
sented  above  are  an  accurate  representation  of  viscosity  and  non-equilibrium 
rate  kinetics  and  most  urgently,  in  light  of  the  recent  measurements  in  ref.  [37] 
(that  confirmed  earlier  speculations  of  anomalous  ionization)  more  realistic  ion¬ 
ization  models  including  microturbulent  effects. 
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5  Conclusions 

We  have  addressed  in  detail  the  problem  of  anomalous  transport  in  the  MPD 
thruster  and  showed  through  numerical  simulation  how  the  inclusion  of  turbu¬ 
lent  effects  in  fluid  codes  can  improve  the  predictability  and  realism  of  fluid 
codes.  We  discussed  the  theoretical  badcground  and  derivations  for  the  trans¬ 
port  models,  and  reduced  our  models  to  the  form  of  polynomials.  We  included 
these  polynmials  in  an  advanced  fluid  code  and  studied  various  features  of  the 
flow  that  are  impacted  by  the  effects  of  plasma  turbulence. 

Now  that  the  models  have  been  derived  and  tested  we  have  proceeded  to 
the  dissemination  of  the  results  to  other  workers  in  the  field.  Since  we  have 
cast  our  bottom-line  results  in  the  form  of  simple  polynomials  that  depend  only 
on  macroscopic  properties,  they  are  thus  ready  for  inclusion  in  any  two-fluid 
code  which  could  then  be  used  for  more  realistic  calculations  of  MPD  thrusters 
including  the  effects  of  plasma  turbulence  caused  by  the  existence  of  instabilities. 
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Figure  17:  Electron  Hall  parameter  in  regions  where  anomalous  transport  is 
active  {uie/vti  >  1.5).  Critical  regions  are  those  with  higher  values  of  the  Hall 
parameter. 


Figure  18;  Ratio  of  “anomalous”  to  “classical”  resistivities  (i.e.  using  the  re¬ 
sistivities  calculated  from  the  anomalous  and  classical  runs).  (J  =  18  kA, 
rh  =  6  g/s). 


Figure  19:  Heavy  species  temperature  for  the  simulation  with  classical  transport 
(J  =  18  kA,  m  =  6  g/s). 
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Figure  20;  Heavy  species  temperature  for  the  simulation  with  anomalous 
transport  showing  regions  with  enhanced  ion  heating  due  to  microturbulence. 
(J  =  18  kA,  rh  =  6  g/s). 


Figiire  21:  Ratio  of  ion  to  electron  temperatures  for  the  simulation  with  classical 
transport  (( J  =  18  kA,  m  ==  6  g/s)).  Note  that  this  ratio  nowhere  exceeds  unity. 


Figure  22:  Ratio  of  ion  to  electron  temperatures  for  the  simulation  with  anoma¬ 
lous  transport  (J  =  18  kA,  m  =  6  g/s). 
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6  Appendix 


The  elements  of  the  dispersion  tensor  Aj  needed  for  the  calculations  carried  in 
this  report  were  derived  in  refs.  [5]  and  are  quoted  here  for  the  sake  of  com¬ 
pleteness. 
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where  /„  is  the  modified  Bessel  function  of  the  first  kind  of  order  n  and  where 
we  have  used  the  following  definitions 

_  w  -f  nuct  -  kiUie  -t-  iut .  ^  ~  ks_uj 

tj  =  W  +  tl/e;  Oo  =  T:^  =  CeO+^^:  (99) 

KfVtg  Kg  Vfg 

—  j^(Cen  ~  0o)»  ^en  =  ^  (Cen)  >  =  fn(Me)-  (f09) 

where  Z  is  the  standard  plasma  dispersion  function  and  where  the  thermal 
velocity,  plasma  frequency  and  cyclotron  frequency  of  species  a  are,  respectively, 
given  by 

e,,  =  ;  u;„  =  (101) 

\  «om,  J  m, 

7  Non-technical  Material 

7.1  Publications  supported  by  the  AFOSR  Grant 

The  following  four  papers  from  the  AFOSR-supported  work,  listed  in  the  refer¬ 
ence  list  and  refered  to  here  with  a  citation  number,  have  already  been  published: 
refs.  [4,  38,  27,  28). 

A  paper  entitled  “Mass  Savings  Domain  of  Plasma  Propulsion”  by  E.Y. 
Choueiri,  A.J.  Kelly  and  R.G.  Jahn  has  been  accepted  in  the  Journal  of  Space¬ 
craft  and  Rockets  and  will  appear  in  the  issue  of  Octoberl993. 

A  paper  entitled  “Numerical  Fluid  Simulation  of  an  MPD  thruster  with  a 
Real  Geometry”  by  G.  Caldo,  E.Y.  Choueiri,  A.J.  Kelly  and  R.G.  Jahn  will 
be  presented  at  the  upcoming  International  Electric  Propulsion  Conference  in 
Seattle,  September  1993  under  paper  number  IEPC-93-072. 

Three  papers  are  currently  being  prepared  for  submission  to  the  following 
three  journals:  Journal  of  Plasma  Physics,  The  Physics  of  Fluids  and  Journal 
of  Propulsion  and  Power. 

The  titles  are,  consecutively:  “Anomalous  Transport  from  a  finite-beta  col- 
lisional  instability”,  “A  finite-beta,  collisional  instability:  Theory  and  experi¬ 
ments”  and  “Plasma  Fluid  Simulations  of  MPD  Thruster  Flows  with  anomalous 
TVansport”.  The  authors  are:  Choueiri,  Kelly  and  Jahn  for  the  first  two  papers 
and  Choueiri,  Caldo,  Kelly  and  Jahn  for  the  last. 

Finally  the  following  eight  technical  reports  on  work  conducted  under  this 
grant  have  been  published  worldwide  by  the  Mecahnical  Aeropsace  Engineering 
Department  of  Princeton  University  between  Feb  1st  1991  and  Feb  1st  1993. 
They  are  listed  in  the  reference  list  and  refered  to  here  with  a  citation  number: 
ref.  [39,  40,  41,  42,  43,  44,  45,  46] 
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7.2  Particpating  Professionals 

•  Dr.  Edgar  Y.  Choueiri,  first  as  a  graduate  student  and  research  assistant 
and  subsequently  as  a  research  associate. 

•  Dr.  Arnold  J.  Kelly,  Senior  Research  Engineer. 

•  Dr.  Robert  G.  Jahn,  Professor 

•  Giuliano  Caldo,  Graduate  Student  and  Researdi  Assistant 

7.3  Advanced  Degrees 

A  substantial  portion  of  the  Ph.D.  thesis  of  Dr.  Choueiri  consisted  of  work 
conducted  under  this  grant.  His  thesis[5]  bears  aknowledgement  to  AFOSR. 
His  degree  was  awrded  in  October  of  1991. 

A  substantial  portion  of  the  Masters  thesis  of  Mr.  Giuliano  Caldo  consisted 
of  work  conducted  under  this  grant.  His  thesis  is  currently  in  preparation. 

7.4  Professional  Interactions 

We  have  presented  the  results  of  this  work  at  three  conferences:  International 
Electric  Propulsion  Conference  1991,  Joint  Propulsion  Conference  1992  and 
Symposium  on  Nuclear  Power  and  Space  1993.  We  will  be  presenting  some 
of  our  final  results  at  the  upcoming  IPEC  conference  in  September  of  this  year. 

7.5  Other  Relevant  Information 

One  of  the  goals  we  set  for  this  work  from  the  outset  is  the  developement  of 
anomalous  transport  models  that  can  be  easily  used  by  other  workers.  Although 
the  methodology  and  derivations  are  quite  complex  we  have  strived  toward  sum¬ 
marizing  the  results  in  the  form  of  compact  mathematical  formulae.  The  re¬ 
sulting  expressions  can  be  used  by  the  numerical  analyst  without  much  detailed 
familiarity  with  the  mathematics  and  physics  that  lay  behind  the  derivations. 
The  detailed  documentation  offered  here  and  in  our  publications  allow  the  more 
interested  reader  to  go  deeper  in  the  theoretical  aspects  of  the  problem. 

We  have  been  transmitting  the  results  of  this  work  to  other  workers  in  the 
field  especially  to  colleagues  at  NASA  Lewis,  NASA  JPL,  Philips  Lab  at  Ed¬ 
wards  Air  Force  Base,  MIT,  Stuttgart  Univresity  and  the  University  of  Pisa. 
Some  researchers  at  some  of  these  institutions  (e.g.  University  of  Pisa)  have 
already  implemented  our  transport  models  in  their  codes  and  others  have  told 
us  about  their  plans  to  use  them  in  their  codes. 

Finally,  we  would  like  to  stress  that  this  work  is  the  first  of  its  kind  in  plasma 
propulsion  and  we  believe  it  has  chartered  a  new  territory  for  the  modelling 
of  plasma  thrusters.  The  techniques  developed  here  can  be  used  to  address 
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another  important  problems  of  MPD  thruster  research,  most  notably  the  role 
of  turbulence  in  ionization. 
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